Unquenched determination of the kaon parameter Bk from improved staggered 

fermions 



Elvira Gamiz, Sara Collins, Christine T.H. Davics 
Department of Physics & Astronomy, University of Glasgow, Glasgow, G12 8QQ, UI^\ 

G. Peter Lepage 

Laboratory of Elementary-Particle Physics, Cornell University, Ithaca, New York 14853, USA 

Junko Shigemitsu 

\Q • Physics Department, The Ohio State University, Columbus, OH 43210, USA 

o . 

Matthew Wingate 

^ , Institute for Nuclear Theory, University of Washington, Seattle, WA 98115, USA 



(N 



HPQCD and UKQCD collaborations 



The use of improved staggered actions (HYP, Asqtad) has been proved to reduce the scaling 
corrections that affected previous calculations of Bk with unimproved (standard) staggered fermions 
in the quenched approximation. This improved behaviour allows us to perform a reliable calculation 
of Bk including quark vacuum polarization effects, using the MILC configurations with nj =2 + 1 
| flavours of sea fermions. We perform such a calculation for a single lattice spacing, a = 0.125 fm, 

, and with kaons made up of degenerate quarks with m s /2. The valence strange quark mass m 3 is 

fixed to its physical value and we use two different values of the light sea quark masses. After a chiral 
£C) ' extrapolation of the results to the physical value of the sea quark masses, we find Bk = 0.83 ± 0.18, 

where the error is dominated by the uncertainty in the lattice to continuum matching at O(ct^). The 
• matching will need to be improved to get the precision needed to make full use of the experimental 

' data on ek to constrain the unitarity triangle. 

PACS numbers: 

<L> 

I. INTRODUCTION 

In the last few years lattice calculations have started achieving the level of precision and the control of uncertainties 



necessary to extract phenomenologically relevant results Q] . Simulations that incorporate quark vacuum polarization 
effects are required for this task, since the uncontrolled errors associated with the quenched approximation are usually 
the main source of uncertainty in these calculations. An example of this can be seen in the study of indirect CP 
violation in the neutral kaon system. 

The CP violating effects in K° — K° mixing are parametrized by Ek, which is defined in terms of the decay 
amplitudes of the physical neutral kaon states Kl and K$ into a state of two pions with isospin equal to as 

_ A(K L -> (7T7r) J=0 ) 
£K A(K s - (7T7r) J=0 ) ■ W 

Experimentally, this quantity is known with a precision of a few percent. On the other hand, the theoretical calculation 
in the Standard Model, under reliable assumptions, yields the expression 

£ ^7fW mMl2 ' (2) 



'Electronic address: e.gamiz@physics.gla.ac.uk,c. davies@physics.gla.ac.uk 



2 



where AMk is the mass difference between Kl and Kg, a very well measured quantity, and Myi is defined by 

2m K M* 12 = (KZ\H* s f 2 \K a ) . (3) 

Ek is thus determined by the hadronic matrix element between K° and K° of the AS — 2 effective Hamiltonian 

H^ff 2 = C AS = 2 C(p) J d 4 xQ AS=2 (x) (4) 

with 

QaS=2(x) = [Sa-ffidalv-A (x) [s b j^d b ] v _ A (x) . (5) 

The Wilson coefficient C(p) is a perturbative quantity known to NLO in a s in both the Naive Dimensional Reg- 
ularization (NDR) and the 't Hooft-Veltman (HV) schemes. The coefficient Cas=2 includes known functions of the 
masses of particles that have been integrated out, and also depends on a certain combination of Cabibbo-Kobayashi- 
Maskawa (CKM) matrix elements. It is these latter elements about which we would like to obtain information. The 
combination of CKM parameters entering in Cas=2 is Im (Vt s V t * d ) , which is equivalent to the unitarity triangle (UT) 
combination of parameters fj [(I — ~p) + const.] |2(. For reviews where explicit expressions for the coefficients in l|4"| 
can be found, see 0. The two-loop expression for the Wilson coefficient C(/x) is given by 



run ji + ^g^4 
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*MsMr /0 ° > (6) 



where 70, 71 are the AS = 2 anomalous dimension at one- loop and two-loop respectively, and /3q, 0i are the first two 
coefficients of the QCD beta function. 70 = 1 is independent of the number of flavours rif, while the other quantities 
depend onnj in the following way 

A-i(u-25£); A -§(„»- 7 ~.i(-r + ^) ; ,r, 

where only the NLO coefficient 71 depends on the choice of scheme and it is given in the MS(NDR) scheme we are 
going to use to quote the results in this work. The other coefficients are universal. Although C(p) depends analytically 
on the number of flavours through the parameters in Q , the numerical dependency on n f in the conversion from the 

WS{NDR) scheme value B K m - NDR (2GeV) to B K is negligible j|. 

The matrix element (K°\Q as=2\K°) , that encodes the non-perturbative physics of the problem, is usually normal- 
ized by its Vacuum Insertion Approximation (VIA) value, defining Bk as the ratio 

B (u)= _J^\^=2(^)\K°) ( ) 

KW |(^|s7o75rf|0)(0|s7075rf|^ )' 

The renormalization group invariant form of Bk, the so-called Bk, is often used to quote results. It is defined as 

B K = C(fi)B K (fx) , (9) 

with C(/i) the Wilson coefficient in the effective Hamiltonian (0J, given at two- loops in ©. 

One can try to constrain the value of the combination of CKM matrix elements in Cas=2 using the experimental 
value of ek and a theoretical calculation of Bk- That constraint gives a hyperbole in the p — 77 plane, where p and 
r\ are the usual UT parameters . The largest source of uncertainty in the final results comes from the error in the 
determination of Bk IHEJ- Improvement in the calculation of this kaon parameter, to reach at least the same level 
of precision (a few percent) as the rest of the errors entering in the analysis of sk, is thus crucial in order to get 
information about the UT. In addition, ek can be a very powerful probe of new physics 0. 

Due to the phenomenological relevance of Bk its calculation has been addressed many times using different tech- 
niques. The three main continuum QCD-based approaches that have been used in its determination are QCD-Hadron 
duality y|, three-point function QCD Sum Rules H and the 1/N C (N c — number of colours) expansion Hflim. Many 
of these determinations give results for Bk only in the chiral limit, since that is a very well established limit in the 
continuum that simplifies the calculation considerably. Chiral corrections however are essential to get a value of Bk 
useful for phenomenologists, since they represent more than 50% of the final number In fact there is only one 
recent determination of the physical Bk pT| , although there is work in progress to calculate chiral corrections to this 
parameter in the framework of the 1/N C expansion fl^.ll3|. 
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The fourth technique that has been used to calculate Bk is lattice gauge theory. It offers model independent results 
and has the potential to reduce the error to the level required by phenomenology. Reviews of lattice determinations 
of Bk in the quenched approximation and some preliminary results including quark vacuum polarization effects can 
be found in Q. 

So far, the value of Bk that the phenomenologists have been using in their studies of the UT Q is a lattice result 
by the JLQCD collaboration [lfij which uses unimproved staggered quarks in the quenched approximation. That is 

the most complete analysis of Bk with lattice techniques to date. The value given in 

M is B K IS - NDR {2 GeV) = 

0.628 ± 0.042, which corresponds to Bk = 0.86 ± 0.06. This number is the result from the extrapolation to the 
continuum limit of the values obtained for seven lattice spacings, with both gauge invariant and non-invariant operators 
-see section III CI and the error includes an estimate of order a 2 and a 2 corrections. Finite volume effects were also 
studied for two different values of a. 

The main source of uncertainty in the JLQCD calculation (not included in the quoted uncertainty) is the unknown 
error from quenching, which could be as large as 15% according to the Chiral Perturbation Theory (CHPT) estimate 
performed in 16]. In order to have a prediction at a few percent level it is thus necessary to perform unquenched 
calculations of Bk that eliminate the quenching uncertainties. Such calculations have now become feasible. 

Another feature of the calculation in [l£j is that it is affected by large scaling violations. We will see that the scaling 
behaviour is much better using improved staggered actions instead of the standard unimproved staggered action used 
by the JLQCD collaboration. 

A third way of improving the JLQCD calculation would be to incorporate SU(3) breaking effects by using kaons 
made up of non-degenerate quarks, instead of degenerate quarks with m s /2. However, these effects arc not expected 
to be important as first estimated using CHPT bySharpe (the effects were estimated to be ~ 5%) and as the 
preliminary results from unquenched domain wall 17] and staggered LIS] quarks seem to indicate (the authors found 
differences ~ 3% between the degenerate and non-degenerate results). 

The goal of this work is to perform a calculation of Bk including quark vacuum polarization effects, which will 
eliminate the irreducible and unknown quenching uncertainty. In this calculation we will use improved staggered 
fermions that have been proved to reduce the large 0(a 2 ) discretization errors generated by the taste-changing 
interactions. 

Preliminary results from this study were presented in 19]. The matching coefficients needed in the calculation 
of the renormalized Bk with the action used in our unquenched simulations were not available at that time, so an 
approximate renormalization was performed in order to get those preliminary results. The correct renormalization 
coefficients have since been calculated |2£j and have been used to obtain the results reported in |2l| and in the present 
article -see next section for more details about the renormalization pro cess. 

Other recent preliminary unquenched results can be found in [jj], u3> US Hj> u3 • 

II. THEORETICAL FRAMEWORK 

In this section we briefly describe the fermion formalism we are going to use in this work, the staggered formalism, 
as well as the fermion and gauge actions used in the quenched and unquenched calculations reported here. We also 
set up the notation for the lattice operators needed in our study and justify our choice of external states. 

A. Staggered fermions 

When the quark action is written in terms of staggered fermions, it becomes spin diagonal and we can drop three 
of the four components of the staggered field x- The standard (unimproved) form of the staggered fermion action is 
thus 



gunim. ^ ^ 



2 



(10) 



where n = (ni, n%, n^) parametrizes the lattice site, rj^{n) = (— l)" lH ^'V-i and x is a 2>{colours) x l(spin) 
component object. The relation between the staggered field x an d the naive fermion field \& at each lattice site is 

*(n)=r nX (n) *(n)=x(n)rt, (11) 

with 

r„ = ( 7 i) ni (72)" 2 (73)" 3 (74)" 4 . (12) 
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The staggered action has a remnant of chiral symmetry which ensures that the Goldstone boson pion mass vanishes 
at ma = 0. But the main advantage of staggered fermions is that they are computationally very efficient and 
unquenched simulations with three flavours of sea quarks are possible at present with lighter sea quark masses than 
those achieved with other fermion formulations. 

The fundamental disadvantage suffered by staggered fermions is the fact that each flavour field comes in four different 
tastes. In the continuum limit the four tastes are degenerate and extra copies can be removed by hand. For sea quarks 
this involves taking the fourth-root of the quark determinant in unquenched simulations and, in perturbation theory, 
the division of each fermion loop by a factor of 4. The validity of the fourth root procedure has not been rigorously 
proven but tests of it are encouraging [2(j. At non-zero lattice spacing the situation is more complicated due to 
the existence of interactions of <D{a 2 ) that violate the taste symmetry and are potentially dangerous. However, the 
negative effects of these taste-changing interactions can be reduced by adding improvement terms to the action as 
pointed out in [27|]. This issue is discussed in more detail in the following section. 



B. Improved staggered actions 



Large Ola 2 ) discretization errors have been found in the calculation of masses and matrix elements using unimproved 
staggered fermions. The origin of these anomalously large discretization errors, as well as other effects, such as the 
large size of perturbative corrections, are the taste changing interactions that break the taste symmetry at nonzero 
lattice spacing. These (unphysical) interactions can be systematically removed using the method of Symanzik and 
replacing the gauge link U in I|1UI) by a fat link V, that is a weighted combination of different staples. For studies of 
the effects of using fat links see |28|, |2{| |3(j . 

In the quenched simulations in Section llVl we use two different improved actions, the HYP and the Asqtad. The 
construction of the HYP [3l| action involves three levels of APE smearing with projection onto SU(3) at each level. 
The smearing is restricted in such a way that each fat link includes contributions only from thin links belonging to 
hypercubes attached to the original link. 

The Asqtad action we analyze here and use in our unquenched simulations has been extensively used in the past for 
simulating light sea and valence quarks It nas been designed to reduce the taste symmetry breaking effects and 
remove all other 0{a 2 ) discretization errors. The virtue of this improved action is that it allows for precise calculations 
with light sea quarks. 

Together with the fat links, the other ingredients that differentiate the fermion part of the Asqtad action from the 
unimproved one are the so-called Lepage corrections and the Naik term. The Asqtad action, written in terms of four 
component naive fermions, has the form 



cAsqtad 4 
S f = a 



^ lll a \ » 6 " 



*(x) 



(13) 



with 



73— link 



tadpole improved 



- { -3 [UUU + 3a,,) - U^U^U^ V(x - 3a p 



\U*(x + a,)-U^(x-a„ 

u 



(14) 



and being the usual covariant derivative with the thin link variable replaced by an updated variable 
V^(x) = Vfi(x) — Sp^p — i~ U^(x). In this expression V l _ l (x) is the fat link 

I symmetrized ^C 2 -) (15) 



and the second term is the Lepage term that removes a low momentum 0(a 2 ) error. For the exact definitions of the 

7 *. (2) 



derivatives V* and W see, for example, |3 
The Asqtad fermion action is coupled in tl 
Symanzik improved after tadpole improvemei 
gluon action, we use the usual unimproved Wilson glue action 



The Asqtad fermion action is coupled in the un quen ched simulations to an SU(3) gluonic action that is one-loop 
Symanzik improved after tadpole improvement |34l l35|. In the quenched calculations, instead of using the improved 
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C. Definition of the operators 



To construct the four-fermion operators, we collect the stag gere d fields x into a set of Dirac fields q(2N) that live 
on the even lattice sites and are spread over a unit hypercube |36fl 

q{2N) ai = ^^2(T A ) a ix(2N + A), (16) 

A 

where a and i are the Dirac and taste indices respectively, and with A running over the vertices of a hypercube 
(A^ = or 1, fj, = 1, . . . ,4). The bilinear quark operators with spin structure 75 = Ts and taste structure £t = 
are defined by 

O ST = q(2N)( ls ® HT)q(2N) = -L £ ^( 2 7V + A) X (27V + B)~tr (r^lsl^) (17) 



The four-quark operators can be build with the bilincars in <|17[1 considering two different contractions of the colour 
indices 



01 = q a h Sl ® 60<Z 6 ■ q b (ls 2 ® fe)g a = ( ^) X! ^I tr ( r A7Sir B r^) X ^ • X^tr (r^r^) Xd , 

^ ' ABCD 

2 = q a h Sl ® 609° • g- fc ( 7 s 2 ® 60<z 6 = (i) E xijtr (r^r^) x % ■ x^jtr (r^r^) x h D (18) 



ABCD 



where we have suppressed the hypercube label 2N for simplicity. The operators in l|18|) are known as one-colour-trace 
and two-colour-trace operators respectively. 

To make these operators gauge invariant we insert gauge link factors connecting the quark fields according to 



^ ' ABCD 

O, = (yX £ Xa^^TbT^)^ -xl^ (19) 



ABCD 



Here, we suppress again the hypercube label 2N for simplicity and add superscripts n f i to label the different continuum 
flavours. Notice that we consider fat links only in the action and those that we introduce in the operators are thin 
links. The only improvement in the operators we consider is tadpole improvement. This we carry out by dividing 
the gluon fields by appropriate factors of the mean link uo defined as the fourth root of the average plaquette. The 
operators used in the calculation with improved actions are thus exactly the same as those used in the unimproved 
simulations. In this work we also use gauge non-invariant operators, that do not incorporate gauge link factors but 
instead are calculated on gluon configurations fixed to Landau gauge. 

In the transcription of the continuum operators involved in the calculation of Bx to the staggered ones we have 
to take into account the taste degree of freedom. We choose the external kaon to have taste structure £t = 75 (and 
thus fix the taste structure of the vector and axial operators with nonzero matrix elements), since only the mesons 
with such structure become massless in the chiral limit. In addition, in order to avoid mixing with operators with 
different taste structure, we follow the two-spin-trace formalism described in references and introduce two sets 

of valence quarks s\, d\ and S2, di for each of the bilinears constituting the four-fermion operator Qas=2- Matrix 
elements are then taken between kaons K® = sidi and K\ = d-is-i- A detailed explanation of the two-spin-trace 
formalism for Bx can be found in [23, US ■ 



III. PERTURB ATIVE RENORMALIZATION 



An important step needed in the determination of Bk (as well as any renormalized quantity) is the calculation of the 
coefficients that match the lattice matrix elements to the continuum ones. The operator Qas=2 in (0 is renormalized 
multiplicatively if one has a regularization with exact chiral symmetry, such as the dimensional regularization in 
the continuum or Ginsparg- Wilson fermions on the lattice. For staggered fermions, although full chiral symmetry 
is broken at non-zero lattice spacing, there is a remnant U(l) symmetry which allows us to have the mixing under 
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control at a given order in perturbation theory. In particular, the one-loop calculation involves only four operators, 
as explained below. 

In this work we use a one-loop perturbative matching to relate the lattice and continuum operators. A general one- 
loop matching of the bare lattice operators O l °" tt ^ to the renormalized continuum operators Of mt ' at an intermediate 
scale [i can be expressed as 



Or'-(M) - 0\ att - + S^l £ (_ 7 « m (H + dj) Of 1 ' + 0(a s (q*f) , (20) 

j 



with "f\ J the one-loop anomalous dimensions matrix and Cy the one-loop matching coefficients, depending on the 
continuum scheme. 

The operator Qas=2 is the product of two V — A currents. The QCD corrections to the bare lattice four-quark 
operators affect the vector and axial parts differently; as a consequence currents of the form V + A are generated. A 
minimal set of lattice operators that matches to the continuum, closes under renormalization and has non-vanishing 
K — K matrix elements is 

V (i) = (s a d b ) v (s b d a ) v , = (s a d a ) v (s b d b )v , 

A« = (s a d b ) A (s b d a ) A , A( 2 ) = (s a d a ) A (s b d b ) A , (21) 

where V and A are the vector and axial currents with taste structure £t = 75- The colour indices a, b indicate 
which fields are connected by gauge link factors, according to l|19|) . The superscripts 1 and 2 indicate that they are 
one-colour-trace or two-colour-trace operators. 

Perturbative calculations and non-perturbative simulations have to be done with exactly the same action and 
operators. The one-loop lattice to continuum matc hing coefficients for the A5 = 2 four-quark operators in l|21|l have 
been calculated for unimproved staggered fermions [3!j , as well as for HYP staggered fermions [4(| ; although this last 
calculation only considered gauge invariant operators. The matching calculation for the Asqtad action, with both a 
Wilson gauge action and an improved gauge action as described in Section ( fllB j) . was done perturbatively at one-loop 
in | 20| for gauge invariant as well as gauge non-invariant operators for a general gauge. 

The calculation in |20| was performed with two independent methods. First, the authors evaluated the corresponding 
one-loop diagrams in two different ways: by directly calculating the various diagrams and by first separating off the 
part which can be inferred from the renormalization of the current operators. Second, the lattice integrals were 
evaluated both algebraically and numerically. For the algebraic evaluation, the authors in 20l | expanded the diagrams 
around the continuum limit. This produces a set of lattice tadpole integrals which one then reduces to a minimal set 
of master integrals using computer algebra. The agreement between these two rather different methods provides a 
strong check on their results. 

The size of the one-loop corrections to the tree level matching for the Asqtad action (with and without improved 
glue) is very similar to what is found with unimproved staggered and other improved staggered actions, such as 
the HYP action. In the calculation of Bk thus the reason to use an improved action is not the reduction of the 
matching factors (which are not large in the unimproved calculation) but the correction of the bad scaling behaviour, 
as we have already pointed out. That is not the case with other weak matrix elements relevant in the study of CP- 
violating effects, for which the use of improved staggered actions greatly reduces the size of the perturbative matching 
coefficients -see JI3 f° r further discussions and references. 

Two comments are in order with respect to the results in [2(j. First, the authors found that the matching coefficients 
obtained for the Asqtad action with improved glue are very similar to those with an unimproved glue action -see 
equations and i)23|) below. They conclude that the improvement in the glue action is not crucial in order to 
reduce the size of the perturbative coefficients, as already suggested in ^lj. Another conclusion from that work is 
that we do not expect anomalously large 0(a 2 s ) corrections using the Asqtad action. 

For completeness, we write here the matrices Cy in J5DJ), taken from 20], that must be used to obtain 
(K°\Q&.s=2(lJ>)\K } MS from the Asqtad bare matrix elements in the Appendix 1X1 They are written in the basis 
(y^,V^,A^\A^). They are 



unimp. 

inv. 



-14.3796 1.3606 -0.0520 -0.0394 

3 -19.8263 -0.0389 0.0263 

-0.0520 -0.0394 -15.1796 3.8005 

-0.0389 0.0263 3 -10.6495 

-12.3228 2.6405 -0.0520 -0.0394 
3 -13.4411 -0.0389 0.0263 



I -0.0520 -0.0394 -12.4428 3.0005 I 
-0.0389 0.0263 3 -12.4811 
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p 


Volume 


Hconfs 


cr x (GeV) (im s /2 


La (fm) 


arriK 


a a (l/a) 


UKQCD n f = 


5.7 


12 3 x 24 


150 


0.837(6) 0.086/0.064 


2.9 


0.417 


0.31 


5.93 


16 3 x 32 


50 


1.59(3) 0.039/0.030 


2.0 


0.220 


0.21 



TABLE I: Parameters in the quenched simulations. The values for m s /2 are for the HYP and Asqtad staggered actions 
respectively. 



for both gauge invariant and gauge non-invariant operators (in the Landau gauge) using an unimproved glue action. 
And 



imp. glue 
inv. 



/-(imp. glue 
noninv. 



-15.1091 
3 

0.5080 
0.4411 

-12.0486 
3 

0.5080 
0.4411 



1.7606 
-19.3177 

0.2406 
-0.1337 

2.9206 
-12.2869 

0.2406 
-0.1337 



0.5080 
0.4411 
-15.8691 
3 

0.05080 
0.4411 
-12.1687 
3 



0.2406 
-0.1337 

4.0806 
-10.8972 

0.2406 
-0.13370 

3.2806 
-11.3269 



(23) 



for gauge invariant operators and gauge non-invariant operators using the improved gluon action. These results are 
obtained after tadpole improvement with the mean link uq defined as the fourth root of the average plaquette. The 
one-loop contribution to this parameter is itg = 7r/3 for the unimproved gluon action and — 0.7671 for the 

improved gluon action, where is defined as uq = 1 — a s UQ + 0(a 2 s ). 

In addition to the matching for the four-fermion matrix element in © we must also account for the renormalization 
of the axial currents in the denominator of ©, which is non- vanishing with the Asqtad action for the definition the 
axial current of l|21|l . The axial current renormalization is multiplicative and can be written in the form 1 + a s /irZA, 
with Z™^ nimp = 1.206, Z™?' = 1.237 when using unimproved and improved glue respectively with gauge invariant 
operators, and Z™™Iimp. — 1-435, Z™£%£ v — 1.291 when using unimproved and improved glue with gauge non- 
invariant operators in the Landau gauge. 



IV. SCALING BEHAVIOUR OF THE IMPROVED STAGGERED ACTIONS 



The first issue we analyze is the impact of using improved staggered actions in the calculation of Bk , in comparison 
with the unimproved staggered action analyzed in |l5j . This study is carried out in the quenched approximation and 
for two different improved actions: the HYP |3l| and the Asqtad with Wilson glue. 

In the next two subsections we describe the results obtained for gauge invariant operators only. We postpone the 
discussion on the differences found between using gauge invariant and gauge non-invariant operators until the last 
subsection IIVCI 



A. Simulation details 



We have two ensembles of 150 and 50 configurations at j3 = 5.70 and {3 = 5.93 respectively, generated using the 
Wilson gluon action for the three staggered fermion actions that we are going to analyze, unimproved, HYP and 
Asqtad. The values of the parameters used in the simulations are shown in Table [I] We choose these parameters to 
be the same as those used by the JLQCD collaboration in order to make a clear comparison with their results. In 
particular, we match kaon masses at a given /3 to those of the JLQCD collaboration (4^, which fixes the strange quark 
mass to the values listed in Table [I] The lattice spacings, determined from m p , are also taken from reference [4*2^ . 
again for consistent c omp arison with that work. For the same reason, we consider kaons made up with two degenerate 
quarks of m s /2 as in |42j . 



Within the conventions and parameters we have described above, the bare values of Bk we obtain are plotted in 
Figure ^ as a function of the timeslice for the three actions and for the two different values of the lattice spacing. 
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(b) 



FIG. 1: Values of the bare gauge invariant Bk at j3 = 5.70 (a) and 5.93 (b) for unimproved (green diamonds), HYP (red 
squares) and Asqtad (blue circles) quarks. Solid lines show results from our fits. 



The /3 = 5.70 and f3 = 5.93 results present a plateau for 3 < t < 10 and 5 < t < 12 respectively, so we make fits to a 
constant over those ranges of values of t. The results from those fits are also plotted in Figure 2] and the numerical 
values summarized in the left hand side of Table [H] Individual values for the four four-fermion bare lattice operators 
in 1)21(1 normalized to the bare lattice value of the numerator in (JHJ are given in Appendix 1X1 





B b £ re (n f = 0) 




B% s - NDR (2GeV) {n/ = 0) 





invariant non-invariant 


13 


invariant non-invariant 




nf — unimproved 




nf = unimproved 


5.7 


0.876(1) 0.895(1) 


5.7 


0.816(1) 0.841(1) 


5.93 


0.734(5) 0.722(3) 


5.93 


0.720(3) 0.747(3) 




nf =0 Asqtad 




nf = Asqtad 


5.7 


0.815(1) 0.824(1) 


5.7 


0.715(1) 0.729(1) 


5.93 


0.699(7) 0.708(4) 


5.93 


0.647(4) 0.673(4) 




n/ = HYP 




ra/ = HYP 


5.7 


0.782(1) 


5.7 


0.648(1) 


5.93 


0.703(5) 


5.93 


0.627(4) 



TABLE II: Bare B K and renormalised B* /s " iVDii (2GeV) for both gauge invariant and non-invariant operators in the quenched 
approximation using unimproved, Asqtad and HYP staggered fermions. The errors quoted in this table are only statistical. 



B. Results 

To convert the lattice results to the MS — NDR scheme we use the one-loop coefficients from [3^, an d 
as appropriate, with the matching scale /i and the scale for a s equal to [i — q* = 1/a. The values for the coupling 
constant a s {q* = 1/a) used in the matching process for the different lattice spacings, are again taken from [l5| to 
provide comparison with that work. They use oijjg_ NDR {l/a) as given by Aj^g — 230MeV. Their values are listed 

m Table HI After this process we obtain the renormalized continuum value B^ S (1/a), that we can run to 2GeV 
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FIG. 2: Scaling of B^ S NDR (2GeV) with a for improved staggered actions compared to the JLQCD unimproved staggered 
results. Results in this Figure are obtained with gauge invariant operators. 

using the two-loop running of the continuum renormalization group 

Bk(mi) = §^|^( M2 ), (24) 

where C(p) is the Wilson coefficient given in with rif = 0. The renormalization group invariant form Bk can be 
similarly calculated using JfjJ) and its definition in ©. 

The results we obtain for B K fS ~ NDR (2GeV) for the different actions and the two lattice spacings are given in the 
right hand side of Table [H] In the unimproved case the numbers for gauge invariant and noninvariant operators 
agree well with those of the JLQCD collaboration for both lattice spacings. The quenched results for gauge 
invariant operators in Table [H] are plotted as a function of the lattice spacing in Figure [21 In this figure a clear 
improvement in the scaling can be seen when using improved actions, in particular in the HYP case. We can quantify 
the improvement in the scaling by extrapolating our results for [3 — 5.70 and (3 = 5.93 to the continuum limit, 
assuming a quadratic dependence on the lattice spacing. In the unimproved case, the result from this exercise, 0.682, 
is incompatible with the JLQCD result 0.628, that incorporates an estimate of 0(a 2 ) and 0(a 2 ) corrections. In fact, 
in the work by the JLQCD collaboration [l5j the authors disregarded the points at /3 = 5.70 and (3 = 5.85 to perform 
the extrapolation to the continuum limit since these points did not exhibit the quadratic dependence in a. However, 

extrapolations of the HYP and Asqtad actions for our two lattice spacings give 0.619 and 0.620 for B^ s ^ NDR (2GeV) 
in the continuum limit, consistent with the JLQCD result in the continuum limit. We expect such improved scaling to 
survive unquenching and this therefore allows us to perform reliable unquenched calculations with only a few values 
of the lattice spacing and even obtain valuable information from simulations at a single scale. 

The reduction of the discretization errors for staggered fermions using improved actions has already been shown 
for other quantities such as hadron masses [2^. In the calculation of Bk it has been recently studied in 0] with the 
HYP action, leading to the same conclusions as in the present article. 

C. Differences between gauge invariant and gauge non-invariant operators results 

The final error quoted by the JLQCD collaboration in its unimproved staggered study of Bk 0] was very much 
enhanced by the differences found between the results obtained with gauge invariant and gauge non-invariant operators. 
These differences have their origin in 0(a 2 s ) and 0(a 2 ) corrections, which were fitted in 0| leading to uncertainties 
of the same size as 3a 2 with a s — a-p^l/a). 

We have studied this issue with the unimproved as well as with the Asqtad action. A similar analysis with the 
HYP action is not possible since renormalization coefficients for gauge non-invariant operators in the Landau gauge 
are not available. The results for the bare Bk as well as the renormalized B K rs ~ NDR (2GeV) obtained with the gauge 
non-invariant operators defined in Section III CI are given in Table ^ Despite the large improvement in the scaling 
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p 


^confs 


Volume 


o _1 (GeV) 


QifThsea. 


arris/2 ay(l/o) 


6.76 


560 


20 3 x 64 


1.605 


0.01/0.05 


0.02 0.47 


6.79 


414 


20 3 x 64 


1.596 


0.02/0.05 


0.02 0.47 



TABLE III: Parameters in the unquenched simulations. The configurations are taken from the MILC collaboration 34]. The 
lattice spacings are taken from [4^ |. Quarks masses are given in the MILC convention, which includes a factor of mo compared 
to the standard convention. 

found for the Asqtad action, the differences between gauge invariant and noninvariant results are of the same size 
as those observed in the unimproved case. That indicates that these differences are dominated not by the 0(a 2 ) 
corrections but by the 0(a 2 ) corrections, as was already pointed out in reference |15| . Since the improvement of the 
action does not lead to a reduction of the perturbative coefficients in the particular case of the calculation of Bk , the 
differences between the results using the two definitions of operators are not reduced by using the improved actions. 
It would be necessary to perform a two-loop matching to reduce the uncertainty associated with the definition of the 
operators. 

Another conclusion from the comparison of the Asqtad results in Table [H] is that those corresponding to gauge 
invariant operators give values more similar to the JLQCD results as a — > 0, which indicates that the 0{a 2 s ) corrections 
are smaller for these operators. We expect thus more accurate results from gauge invariant operators, which is what 
we use in our unquenched calculation, than from the noninvariant ones. 

V. UNQUENCHED VALUE OF B K 

We now incorporate quark vacuum polarization effects in the calculation of Bk using one of the improved staggered 
actions analyzed in the quenched approximation, Asqtad, since the final goal is to eliminate the irreducible systematic 
error associated with quenching that dominates the total uncertainty in previous determinations of Bk- We use the 
Asqtad action because there are configurations for this action generated with sea masses sufficiently small to perform 
a realistic chiral extrapolation to the physical point |45l l4q . In addition, unquenched simulations using this action 
have been successful in describing a wide range of experimental observables with systematic errors of 3% or less Q. 

We performed an unquenched calculation of Bk with the Asqtad action described in Section Hi Bl using the config- 
urations from the MILC collaboration with nj = 2 + 1 sea flavours |3J| . The results reported here correspond to the 
analysis at one lattice spacing with a = 0.125 fm and two different values of the light sea quark masses. The parame- 
ters used in the unquenched simulations are collected in Table iLUl As in the quenched simulations we use degenerate 
quark kaons with the strange quark mass fixed to its physical value, so no extrapolation of the valence quark mass is 
necessary. A difference with the quenched analysis described in the last section is that for the unquenched case we 
only consider gauge invariant operators. 

A first promising sign in our calculation is that there is little contamination from excited states, as can be seen in 
Figure |3 where we have plotted the bare values of Bk as a function of the timeslice. In particular, for (3 — 6.76, for 
which we have better statistics, we obtain an excellent plateau. To get the bare values of Bk given in the Appendix 
lAl we perform fits to a constant over 5 < t < 27 for [3 — 6.76 and over 8 < t < 30 for (3 = 6.79. Results from the fits 
over the same ranges for the different four-fermion bare lattice operators involved in the calculation, normalized to 
the bare lattice value of the numerator in ©, are given in the Appendix lAl 

The conversion of the values of the bare lattice operators to a value for B^! s ~ NDR (2GeV) has been done pertur- 
batively using the 0(a s ) lattice to continuum matching coefficients from [2Cj collected in the matrix (|23|l . In the 
matching process it is most natural to take a s in the V scheme with values for Nf = 3 from the recent 4-loops 
lattice determination in -see Tabic ILTTI The scale for a s is not determined here so we consider various reasonable 
possibilities. One can also optimise the scale /x in equation 1)2011 |4Sj . but we have used pL — 1/a throughout this work. 

The results we obtain for B^ s ^ NDR (2GeV) performing the renormalization with ay(l/a) are shown in Figure^] 
We plot the result as a function of the light sea quark mass and a decrease of the value of Bk with the reduction of 
the sea quark mass can be appreciated in this figure. For the degenerate valence quark case we are analyzing, the 
chiral behaviour of Bk with the sea quark masses is linear, with the same coefficient for the strange, up and down 
sea masses We can thus extrapolate our results to the physical s and u(d) masses, which yields the result 

B K m - NDR (2GeV) = 0.618(18)(19)(30)(130) . (25) 
The first error in (|25fl is statistical, the second is from the extrapolation to the physical values of the sea quark masses, 
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FIG. 3: Figure showing bare gauge invariant values of Bk as a function of the timeslice corresponding to unquenched simulations 
with m sea = 0.01/0.05 (a) and m Bea = 0.02/0.05 (b). 



TTlsea 


glare gMS- NDR ( 2Ge V) 


rif = 2 + 1 Asqtad 


0.01/0.05 
0.02/0.05 


0.785(11) 0.655(9) 
0.815(18) 0.680(13) 



TABLE IV: Values of the bare Bk and the one-loop renormalized B K ~ (2GeV) for the two values of the light sea quark 
masses and the corresponding statistical errors. The first number in the row labelled m s is the mass of the degenerate up and 
down sea quarks and the second one is the mass of the strange sea quark. 



the third one is from discretizations errors and the final one is from the perturbative conversion to the MS — NDR 
scheme. The value in (|25H is equivalent to Bk = 0.83 ± 0.18, with Bk defined in ©. Note that this value of Bk 
is very similar to the previous quenched staggered result in (Bk = 0.86 ± 0.06), so any final conclusion about 
the enhancement or decrease due to the inclusion of quark vacuum polarization effects in its calculation would need 
a significant reduction of the error quoted in (|2"5)l . 

If instead of using a s in the V scheme to perform the renormalization at one-loop, we use as in the MS scheme at 
a scale 1/a, the result we find is 

Bf" A ' M (2GcV) = 0.637(19) (20) (31) (82) . (26) 

This is the same as the result we obtain when using ay with a scale close to 2/a. 

The total errors in i125[) and i|26|l are dominated by the uncertainty associated with possible 0(a 2 ) corrections in the 
lattice to continuum matching process. On the finer MILC ensembles on which we are planning to redo our calculation 
as the next step in the improvement of the work presented here, the perturbative error bar would be reduced since 
the value of a s (l/a) would be smaller. Using the existing MILC configurations with a = 0.093 fm, the perturbative 
error could be reduced from 20% to 14% (or from 14% to 9% if we use a-jjg(l/a) in the matching). In view of 
the improved scaling behaviour we obtained in the quenched approximation within the Asqtad action described in 
Section IIVI results for two lattice spacings will be enough to perform a reliable continuum extrapolation, reducing 
the discretization errors (third error in l|25|) ). However, what is really needed to reduce the final error in l|25|l to a 
few percent level is a 2-loop matching or a non-perturbative matching method, that eliminates the uncertainty on the 
possible 0(a 2 ) corrections. 
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Gauge invariant B K (2GeV): dynamical vs. quenched 

l i i l i 

O Dynamical Asqtad 
— quenched unimproved (jlqcd) 



» 



0.1 0.2 0.3 0.4 0.5 



FIG. 4: Unquenched value of {ZGeV) as a function of the ratio between the light sea quark mass and the (real) 

strange quark mass. The lines represent the quenched results from Hal . Errors on the points are statistical only. 

A. Discussion of the perturbative error 

The naive error associated with the perturbative matching that we quote in (|25l) is just the result of multiplying our 

central value for B^ s ~ NDR (2GeV) by (ay(l/a)) 2 . In the same way, the perturbative error in (|26|) is the product of 

B K IS ~ NDR (2GeV) in that equation and (a-jg^l / 'a)f \ At the lattice spacing we are working, the unquenched values 

of ay are large -see Table Hilt and translate into a ~ 20% error in the result for B% s ~ NDR (2GeV). 

In principle, we do not expect two-loop corrections as large as 20%, since for the Asqtad action there are perturbative 
corrections for both numerator and denominator in (jSJ) that tend to cancel at any order in the expansion. In fact, 
the shift in the value of Bk due to the one-loop renormalization is only 18%. Further evidence in favour of smaller 
0(a 2 s ) corrections is the fact that the difference between the results we obtain doing the matching with a s at q* = 1/a 
and q* — 2/a, that can be taken as an estimate of the uncertainty from the truncation of the perturbative series, 
is around 4%. The perturbative uncertainty estimated by squaring ajjg{l/a), 13%, then seems more realistic. On 
the other hand, we believe that taking the two-loop error to be just the square of the one-loop shift could be an 
underestimate of the error in the absence of other information. 



VI. SUMMARY AND CONCLUSIONS 

Most of the previous lattice calculations of Bk, in particular those used in the UT analysis, were performed in the 
quenched approximation. This induces a large, essentially unknown and irreducible systematic error into the result. 
Precise simulations with sea quarks are necessary in order to be able to make full use of the experimental data on 
Sk to constrain the CKM matrix. These unquenched simulations are feasible with present computers using staggered 
fermions at light sea quark masses. However, the unimproved staggered action suffers from large taste-changing 
interactions that generate important scaling corrections, as those found by the JLQCD results. We have shown in 
Section llVI that the scaling behaviour is much better when using improved staggered actions. 

This reduction of the discretization errors, together with the existence of unquenched configurations with relatively 
small sea quark masses, makes improved staggered actions an ideal choice for accurate calculations of Bk that 
incorporate light quark vacuum polarization effects. As a first step in this study we have calculated Bk with the 
Asqtad action in two ensembles at a = 0.125 fm and with two different light sea quark masses. In doing that, we have 
used the recent results for the one-loop matching coefficients in J2jj. We obtain, using ay (1/a), 

B~MS- NDR (2GeV) = 0.618 ± 0.136 or, equivalently, B K = 0.83 ± 0.18 , (27) 
or, performing the perturbative matching with otjjg(l/a), 

B^- NDR {2GeV) = 0.637 ± 0.092 or, equivalently, B K = 0.85 ± 0.12 . (28) 
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Both results are compatible within errors with other preliminary unquenched determinations [T?l Il8l I22I 123I |24| . 

The error needs to be reduced further. Some reduction can be achieved by working on finer lattices and incorporating 
staggered chiral perturbation theory results [42|, but a large reduction will require matchings coefficients calculated 
beyond 0(a s ). 

Another issues we would like to investigate in the future are the impact of SU(3) breaking effects and the chiral 
limit value of this quantity that could be compared to recent continuum calculations [lfT |. 
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APPENDIX A: DETAILS OF CALCULATION RESULTS 



In this Appendix we collect the different four-fermion bare lattice matrix elements involved in the calculation of 
Bk for all the cases we analyze in this work. 



£A 



(i) 



A 



(i) 



E A 



(2) 



E^ 



(i) 



u 



(i) 



(2) 



C2) 



Unimproved nf =0 (gauge invariant operators) 



5.70 
5.93 



0.026(1) 0.226(1) 0.012(1) 0.775(1) -0.143(1) 
0.217(5) 0.332(2) 0.091(2) 0.809(2) -0.551(7) 



-0.010(1) -0.009(1) -0.002(1) 
-0.079(2) -0.065(1) -0.0182(4) 



5.70 
5.93 



Unimproved n/ = (gauge non-invariant operators) 



0.035(1) 0.276(1) 0.014(1) 0.772(1) -0.176(1) -0.013(1) -0.010(1) -0.002(1) 
0.211(3) 0.346(1) 0.080(1) 0.800(2) -0.571(5) -0.077(1) -0.056(1) -0.015(1) 



TABLE V: Values of the gauge invariant and gauge non-invariant bare matrix elements in the quenched approximation, shown 
as a ratio with the denominator in for unimproved staggered fermions. 



/Hea (1) 4 1} EA (2) 4 2) E^ t (1) vt E^ (2) vj 2) 

HYP nf =0 (gauge invariant operators) 



5.70 
5.93 



0.048(1) 0.262(1) 0.016(1) 0.772(1) -0.285(2) -0.0194(4) -0.0079(3) -0.0016(1) 
0.181(8) 0.333(4) 0.065(6) 0.798(5) -0.540(13) -0.061(3) -0.042(3) -0.009(1) 

HYP nf — (gauge non-invariant operators) 



5.70 
5.93 



0.041(1) 0.279(1) 0.014(1) 0.769(1) -0.278(1) -0.017(0) -0.008(0) -0.001(0) 
0.164(8) 0.344(4) 0.053(5) 0.796(5) -0.550(1) -0.058(2) -0.039(3) -0.008(7) 



TABLE VI: Values of the gauge invariant and gauge non-invariant bare matrix elements in the quenched approximation, shown 
as a ratio with the denominator in for HYP staggered fermions. 
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\EA 



(i) 



A 



(i) 



(2) 



A 



(i) 



(i) 



(2) 



^4 



(2) 



Asqtad gauge invariant operators 



5.7 

5.93 







0.049(6) 0.285(1) 0.018(1) 0.781(1) -0.278(1) -0.020(1) -0.016(1) -0.004(1) 
0.198(6) 0.351(3) 0.077(4) 0.808(4) -0.596(11) -0.074(2) -0.053(2) -0.012(1) 



n f = 2 + 1 msea = 0.01/0.05 m val = 0.02 



6.76 0.248(6) 0.316(5) 0.076 (2) 0.885(16) -0.605(11) -0.093(2) -0.035(1) -0.0092(4) 



n/ = 2+ l m aea = 0.02/0.05 nw = 0.02 



6.79 0.244(8) 0.323(8) 0.080(3) 0.940(26) -0.630(17) -0.096(2) -0.036(1) -0.0099(4) 



Asqtad gauge non-invariant operators 



n/ 







5.7 
5.93 



0.047(1) 0.284(0) 0.017(0) 0.777(1) -0.267(1) 
0.164(4) 0.338(2) 0.066(2) 0.797(3) -0.548(7) 



-0.018(0) -0.011(0) -0.002(0) 
-0.061(2) -0.041(2) -0.009(4) 



TABLE VII: Values of the gauge invariant and gauge non-invariant bare matrix elements in the quenched approximation and 
with n/ = 2 + 1 sea flavours, shown as a ratio with the denominator in JHJ for Asqtad staggered fermions. 
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